Supplementary material for Warming at the population level: effects on age structure, density, and generation cycles.

library(ggplot2)
library(pander)
library(WaveletComp)

setwd("~/Google Drive/Plodia_populations/Analysis")  #NB set your working directory to whatever is appropriate

adults<-read.csv("Dead adult counts.csv")

adults$temp<-factor(adults$temp)

larvae<-read.csv("Larval counts.csv")

larvae$temp<-factor(larvae$temp)
larvae$food<-relevel(larvae$food, ref="good")
larvae$X5th<-larvae$X5F+larvae$X5M  ## Combine male and female 5ths to give total
larvae$Total<-larvae$X3rd+larvae$X4th+larvae$X5th  ## A few Totals don't equal the sum of 3rds, 4ths, 5ths


larvae2<-subset(larvae, number>10)

Design

• Three temperatures: 27, 30 and 33 degrees
• Two food treatments: good food and poor food (the latter cut with indigestible methylcellulose)
• Three populations maintained at each treatment combination so 18 in total.
• Adult population data are weekly counts of dead adults
• Larval population data are counts of 3rd, 4th and 5th instar larvae from a 3-week old section of food, taken once a week.
• Adult data runs for 82 weeks
• Larval data are available for the first 60 weeks or so

Persistence

NB this analysis and the remaining analyses of population data are carried out using truncated time series with the first 10 weeks removed so as to avoid bias from transients associated with the initial setup of the populations.

Figure S1 shows time series for each population. All populations at 27 and 30 degrees persisted for the whole 82 weeks, as did two of the 33 degree populations, both on poor food. The populations that went extinct did so in less than 6 months, so the high temperature populations either went extinct quite rapidly or persisted for a long time.

rep<-ifelse(adults$box<7,1,ifelse(adults$box>12,3,2))

adults$food<-relevel(adults$food, ref="good")


adults<-data.frame(adults,rep)
rm(rep)
adults2<-subset(adults, week>10 & week<82)

temp_labels <-c("27" = "27°", "30" = "30°", "33" = "33°")
food_labels <-c("good" = "Standard diet", "bad" = "Poor diet")
rep_labels <-c("1" = "Replicate 1", "2" = "Replicate 2", "3" = "Replicate 3")

popplot<-ggplot(data=adults2, aes(x=week, y=total))
popplot<-popplot + geom_line(size = 0.2) + xlab("Time (weeks)") + ylab ("Total adults")
popplot<-popplot+facet_grid(rep~temp+food, 
              labeller = labeller(rep = rep_labels, temp = temp_labels, food = 
              food_labels))+theme_bw()+theme(panel.grid.major=element_blank(), 
              panel.grid.minor=element_blank(), 
              strip.background =element_rect(fill="white", colour = "grey80"))
popplot

Figure S1 Time series for adult counts for each population. The first 10 weeks of data are not shown.

Population statistics

Adults

Figure S2 shows the mean population sizes (as weekly counts of dead adults) for each population. There are very clear differences between the temperatures, with the 30 degree treatments having the highest numbers of adult moths. The effect of food quality appears to be much smaller, if there is an effect at all. Of the 33 degree treatment populations, the two poor-quality food populations which persisted for the duration of the experiment can be seen to have much lower population sizes than the populations maintained at lower temperatures. The one good-quality population which has a mean population size which is similar to the two which persisted is a population which lasted longer than the other high-temperature ones but which eventually went extinct at around 26 weeks.

attach(adults2)
Mean.sizes<-tapply(total,box,function(x) mean(x, na.rm=TRUE))
Temperature<-factor(temp[1:18])
Food<-factor(food[1:18])

Sizes<-data.frame(Temperature, Food, Mean.sizes)

persisted<-c(1,2,3,4,6,7,8,9,10,12,13,14,15,16)

Sizes2<-Sizes[persisted,]

sizeplot <- ggplot(Sizes2, aes(x=Temperature, y=Mean.sizes, fill=Food)) + ylim(0, 165)
sizeplot <- sizeplot + geom_dotplot(binaxis="y", stackdir="center", 
                                    position=position_dodge(0.5)) 
sizeplot <- sizeplot + theme_bw() + scale_fill_brewer(palette="Dark2", 
                                    name = "Food\nquality",
                                    breaks = c("good", "bad"), 
                                    labels = c("Standard", "Poor"))

sizeplot<-sizeplot + ylab("Mean number of adults per week") + xlab ("Temperature °C")
sizeplot

rm(Mean.sizes, Temperature, Food)
detach(adults2)

Figure S2 Mean population sizes for the Plodia populations by temperature and food quality. The first 10 data points for each population have been removed to avoid transients.

We can fit a model to these data to assess effect sizes and the statistical significance of the effects of temperature and food quality. Because most of the high-temperature populations did not persist we’ll restrict the analysis to the 27 and 30 degree treatments only.

Full model

model1<-lm(Mean.sizes~Temperature*Food, 
           data=Sizes, subset=Temperature=="27" | Temperature=="30")

anova(model1)

Analysis of Variance Table

Response: Mean.sizes
                 Df Sum Sq Mean Sq F value   Pr(>F)    
Temperature       1 7851.2  7851.2 42.5928 0.000183 ***
Food              1    6.3     6.3  0.0344 0.857552    
Temperature:Food  1  128.4   128.4  0.6966 0.428150    
Residuals         8 1474.7   184.3                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

rm(model1)

Minimal model

Model selection by removal of non-signficant terms leaves us with only temperature as an explanatory variable

model2<-lm(Mean.sizes~Temperature, 
           data=Sizes, subset=Temperature=="27" | Temperature=="30")

anova(model2)

Analysis of Variance Table

Response: Mean.sizes
            Df Sum Sq Mean Sq F value    Pr(>F)    
Temperature  1 7851.2  7851.2  48.784 3.786e-05 ***
Residuals   10 1609.4   160.9                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

rm(model2)

Since food quality has no significant effect we have a few more df to play with so we can include the two persistent populations at 33 degrees as well.

extinct<-c(5,11,17,18)

model3<-lm(Mean.sizes[-extinct]~Temperature[-extinct], data=Sizes)

summary(model3)

Call:
lm(formula = Mean.sizes[-extinct] ~ Temperature[-extinct], data = Sizes)

Residuals:
    Min      1Q  Median      3Q     Max 
-29.127  -6.317   3.265   6.449  13.183 

Coefficients:
                        Estimate Std. Error t value Pr(>|t|)    
(Intercept)               93.322      4.963  18.802 1.04e-09 ***
Temperature[-extinct]30   51.157      7.019   7.288 1.57e-05 ***
Temperature[-extinct]33  -64.751      9.927  -6.523 4.29e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 12.16 on 11 degrees of freedom
Multiple R-squared:  0.9307,    Adjusted R-squared:  0.9181 
F-statistic: 73.87 on 2 and 11 DF,  p-value: 4.206e-07

rm(model3)

Overall, temperature has a substantial effect on the mean population sizes. Increasing the temperature from 27 to 30 degrees increases mean population size roughly 1.5 times, but going from 27 degrees to 33 degrees means that those populations manage to persist at the higher temperature have a mean size which is roughly a third of the 27 degree populations and a fifth of the 30 degree ones.

Variability of adult populations

attach(adults2)
sds<-tapply(total,box,function(x) sd(x, na.rm=TRUE))

Covars<-data.frame(Sizes, sds)
Covars$CV<-Covars$sds/Covars$Mean.sizes

persisted<-c(1,2,3,4,6,7,8,9,10,12,13,14,15,16)

Covars<-Covars[persisted,]

sizeplot <- ggplot(Covars, aes(x=Temperature, y=CV, fill=Food))
sizeplot <- sizeplot + geom_dotplot(binaxis="y", stackdir="center", 
            position=position_dodge(0.5)) + theme_bw() + scale_fill_brewer(palette="Dark2",
            name = "Food\nquality", breaks = c("good", "bad"), 
            labels = c("Standard", "Poor"))
sizeplot<-sizeplot + xlab("Temperature (°C)")
sizeplot <- sizeplot + ylab("Coefficient of variation for each adult population")
sizeplot

detach(adults2)

Figure S3 Coefficients of variation for each population

30º Populations were slightly less variable than 27º populations. The 33º populations that persisted had much higher CVs but this is possibly more of an effect of their low mean population size than anything else.

Can we discern any statistically significant patterns in the variance of the adult populations?

model1<-lm(CV~Temperature*Food, data=Covars, subset=Temperature=="27" | Temperature=="30")
drop1(model1, test="F")

Single term deletions

Model:
CV ~ Temperature * Food
                 Df Sum of Sq       RSS     AIC F value Pr(>F)
<none>                        0.0099351 -77.159               
Temperature:Food  1 0.0014022 0.0113373 -77.575  1.1291  0.319

model2<-lm(CV~Temperature+Food, data=Covars, subset=Temperature=="27" | Temperature=="30")
drop1(model2, test="F")

Single term deletions

Model:
CV ~ Temperature + Food
            Df Sum of Sq      RSS     AIC F value  Pr(>F)  
<none>                   0.011337 -77.575                  
Temperature  1 0.0090049 0.020342 -72.560  7.1485 0.02547 *
Food         1 0.0011779 0.012515 -78.389  0.9350 0.35882  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

No signficant interaction term and no effect of food. We can include the 33º populations in a final model.

model3<-lm(CV~Temperature, data=Covars)
drop1(model3, test="F")

Single term deletions

Model:
CV ~ Temperature
            Df Sum of Sq      RSS     AIC F value    Pr(>F)    
<none>                   0.013291 -91.436                      
Temperature  2   0.19094 0.204228 -57.186  79.014 2.978e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(model3)

Call:
lm(formula = CV ~ Temperature, data = Covars)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.049353 -0.021899 -0.007172  0.018844  0.056647 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)    0.53089    0.01419   37.41 6.01e-13 ***
Temperature30 -0.05479    0.02007   -2.73   0.0196 *  
Temperature33  0.29838    0.02838   10.51 4.47e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.03476 on 11 degrees of freedom
Multiple R-squared:  0.9349,    Adjusted R-squared:  0.9231 
F-statistic: 79.01 on 2 and 11 DF,  p-value: 2.978e-07

CVs are significantly lower in the 30º populations and then significantly raised in the 33º populations. There is no discernable effect of food quality.

Larvae

Figure S4 shows the mean counts per week for the total number of larvae per food section. By contrast with the data for adults there were fewer larvae in the 30 degree treatment than in the 27 degree one. The two 33 degree treatments which persisted havd very low numbers.

attach(larvae2)
Mean.sizes<-tapply(Total,box,function(x) mean(x, na.rm=TRUE))
Temperature<-factor(temp[1:18])
Food<-factor(food[1:18])

detach(larvae2)

Sizes<-data.frame(Temperature, Food, Mean.sizes)
Sizes<-Sizes[persisted,]

sizeplot <- ggplot(Sizes, aes(x=Temperature, y=Mean.sizes, fill=Food))
sizeplot <- sizeplot + geom_dotplot(binaxis="y", stackdir="center",
          position=position_dodge(0.5)) + theme_bw() + scale_fill_brewer(palette="Dark2", 
          name = "Food\nquality", breaks = c("good", "bad"),
          labels = c("Standard", "Poor"))
sizeplot <- sizeplot + ylab("Mean number of larvae per week") 
sizeplot <- sizeplot + xlab("Temperature °C") + ylim(c(0,17))
sizeplot

Figure S4 Mean number of larvae counted per week per population - compare with Figure S2

Analysis of these data is similar to that for the adult populations. An initial model fitted to the 27 and 30 degree data only is reduced to one with only temperature as an explanatory variable and we can add in the two persisting 33 degree populations to that to give a final minimal adequate model with temperature as an explanatory variable but not food quality. By contrast with the adult populations, larval populations are smaller in the 30 degree treatment than in the 27 degree one, but as with the adult populations the two 33 degree populations which persisted have very low larval counts.

mod1<-lm(Mean.sizes[-c(6,12)]~Temperature[-c(6,12)]*Food[-c(6,12)], data=Sizes)
anova(mod1)
mod2<-lm(Mean.sizes[-c(6,12)]~Temperature[-c(6,12)]+Food[-c(6,12)], data=Sizes)
anova(mod2)
mod3<-lm(Mean.sizes[-c(6,12)]~Temperature[-c(6,12)], data=Sizes)
anova(mod3)

Minimal model

mod4<-lm(Mean.sizes~Temperature, data=Sizes)
anova(mod4)

Analysis of Variance Table

Response: Mean.sizes
            Df  Sum Sq Mean Sq F value    Pr(>F)    
Temperature  2 182.679  91.340  48.464 3.511e-06 ***
Residuals   11  20.732   1.885                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(mod4)

Call:
lm(formula = Mean.sizes ~ Temperature, data = Sizes)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.73718 -0.79327 -0.05288  0.46554  3.14744 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)    13.2564     0.5605  23.653 8.78e-11 ***
Temperature30  -2.6955     0.7926  -3.401  0.00592 ** 
Temperature33 -11.0353     1.1209  -9.845 8.64e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.373 on 11 degrees of freedom
Multiple R-squared:  0.8981,    Adjusted R-squared:  0.8795 
F-statistic: 48.46 on 2 and 11 DF,  p-value: 3.511e-06

rm(Mean.sizes, Temperature, Food, persisted, mod1, mod2, mod3, mod4, Sizes)

Age structure of the larval population

Time averaged proportions of 3rd, 4th and 5th instar larvae are shown in Figure S4. There are clear differences between treatments, with the highest proportions of 5th instar larvae in the 27 degree good food populations and the lowest proportions being in the 30 degree populations.

prop.3rd<-larvae2$X3rd/larvae2$Total
prop.4th<-larvae2$X4th/larvae2$Total
prop.5th<-larvae2$X5th/larvae2$Total

#Calculate mean proportions per box
props.3rds<-tapply(prop.3rd,larvae2$box,function(x) mean(x, na.rm=TRUE))  
props.4ths<-tapply(prop.4th,larvae2$box,function(x) mean(x, na.rm=TRUE))
props.5ths<-tapply(prop.5th,larvae2$box,function(x) mean(x, na.rm=TRUE))


ci.3rds.prop <- tapply(prop.3rd, 
    larvae2$box, function(x) qt(0.975, 
                df = length(x) -1) * sd(x, na.rm = TRUE) / sqrt(length(x)))   
#calculate 95%CI for each mean proportion per box
ci.4ths.prop <- tapply(prop.4th,
    larvae2$box, function(x) qt(0.975, 
                df = length(x) -1) * sd(x, na.rm = TRUE) / sqrt(length(x)))   
ci.5ths.prop <- tapply(prop.5th, 
    larvae2$box, function(x) qt(0.975, 
                df = length(x) -1) * sd(x, na.rm = TRUE) / sqrt(length(x)))   


ordered<-c(1,7,13,2,8,14,3,9,15,4,10,16,6,12)

props<-c(props.3rds[ordered],props.4ths[ordered],props.5ths[ordered])
cis <- c(ci.3rds.prop[ordered], ci.4ths.prop[ordered], ci.5ths.prop[ordered])
Instar<-factor(rep(c("3rd","4th","5th"),each=14))
Population<-factor(rep(c(1,2,3,1,2,3,1,2,3,1,2,3,1,2), times=3))
Food<-factor(rep(c("Standard diet","Standard diet","Standard diet",
       "Poor diet","Poor diet","Poor diet","Standard diet","Standard diet",
       "Standard diet","Poor diet","Poor diet","Poor diet","Poor diet",
       "Poor diet"), times=3))
Food<-relevel(Food, ref = "Standard diet")
Temperature<-factor(rep(c(rep(c("27°","30°"), each=6),"33°","33°"), times=3))

Proportions<-data.frame(props,cis,Instar,Population,Food,Temperature)

## Comment out this line for all instars from 3rd to 5th
Proportions <- Proportions[Instar != "4th",]  

p1<- ggplot(data=Proportions, aes(x=Population, y=props, fill=Instar))

## Remove position = "dodge" for a stacked plot
p1<- p1+geom_bar(stat="identity", position = "dodge") +ylim(0, 0.6)  
#p1 <- p1 + geom_errorbar(aes(ymin = props - cis, 
      #ymax = props + cis), width = 0.1, position = position_dodge(0.9))  

#Comment this out for a stacked plot
p1<- p1+facet_grid(~Temperature+Food,scales = "free_x", 
        space="free_x")+theme_bw() + ylab("Proportion of total larvae") 
p1 <- p1 + xlab("Replicate") + theme(strip.background=element_rect( fill="white"))

#p1<- p1+scale_y_continuous(expand = c(0,0))+
#p1<- p1 + scale_fill_discrete(guide = guide_legend(reverse=TRUE))+scale_fill_brewer()
# test_palette <- c("chocolate3", "steelblue")
test_palette <- c("#ef8a62", "#67a9cf")
# p1<- p1 + scale_fill_brewer(guide = guide_legend(reverse=TRUE))
p1<- p1 + scale_fill_manual(values = test_palette)
p1

Figure S5 Proportions of the total larvae counted in each population that were in the 3rd or 5th instar. 4th instar proportions not shown for clarity.

These differences in the proportions are statistically significant. They can be analysed using a linear model with the mean number in one stadium (either 5th or 3rd instar larvae, it makes little difference) as a response variable and with the total number of larvae as a covariate.

mean.3rds<-tapply(larvae2$X3rd,larvae2$box,function(x) mean(x, na.rm=TRUE))
mean.4ths<-tapply(larvae2$X4th,larvae2$box,function(x) mean(x, na.rm=TRUE))
mean.5ths<-tapply(larvae2$X5th,larvae2$box,function(x) mean(x, na.rm=TRUE))
mean.all<-tapply(larvae2$Total,larvae2$box,function(x) mean(x, na.rm=TRUE))

keep<-c(1,2,3,4,7,8,9,10,13,14,15,16)

Vths<-mean.5ths[keep]  
#Making variables with short names to keep the output readable
Total<-mean.all[keep]  
#Subscripting with "keep" means that only the 27 and 30 degree boxes 
#are being included - because there are only 2 33 degree ones and 
#they are both poor food we can't test for an interaction if these are included.

Food<-larvae2$food[keep]
Temp<-as.factor(larvae2$temp[keep])

mod1<-lm(Vths~Total+Food*Temp)
anova(mod1)

Analysis of Variance Table

Response: Vths
          Df  Sum Sq Mean Sq  F value    Pr(>F)    
Total      1 17.9399 17.9399 139.4633 7.079e-06 ***
Food       1  1.0627  1.0627   8.2611   0.02385 *  
Temp       1  1.0342  1.0342   8.0401   0.02521 *  
Food:Temp  1  0.8885  0.8885   6.9069   0.03401 *  
Residuals  7  0.9004  0.1286                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(mod1)

Call:
lm(formula = Vths ~ Total + Food * Temp)

Residuals:
    Min      1Q  Median      3Q     Max 
-0.3994 -0.2969  0.0303  0.2409  0.3688 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)   
(Intercept)      0.9393     1.1555   0.813  0.44302   
Total            0.4384     0.0831   5.275  0.00115 **
Foodbad         -1.2829     0.3012  -4.260  0.00375 **
Temp30          -1.4403     0.3754  -3.836  0.00640 **
Foodbad:Temp30   1.0899     0.4147   2.628  0.03401 * 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.3587 on 7 degrees of freedom
Multiple R-squared:  0.9587,    Adjusted R-squared:  0.9352 
F-statistic: 40.67 on 4 and 7 DF,  p-value: 6.212e-05

rm(Vths,Total,Food,Temp)

Alternatively, an analysis using the proportion of 5th instar larvae as the response:

propmod<-lm(props~Food*Temperature, data=Proportions, subset=Instar=="5th")
anova(propmod)

Analysis of Variance Table

Response: props
                 Df    Sum Sq   Mean Sq F value   Pr(>F)   
Food              1 0.0153438 0.0153438 10.0619 0.011329 * 
Temperature       2 0.0270404 0.0135202  8.8661 0.007455 **
Food:Temperature  1 0.0010241 0.0010241  0.6716 0.433659   
Residuals         9 0.0137244 0.0015249                    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

propmod2<-update(propmod,~.-Food:Temperature)
anova(propmod2)

Analysis of Variance Table

Response: props
            Df   Sum Sq   Mean Sq F value   Pr(>F)   
Food         1 0.015344 0.0153438 10.4036 0.009091 **
Temperature  2 0.027040 0.0135202  9.1671 0.005476 **
Residuals   10 0.014749 0.0014749                    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(propmod2)

Call:
lm(formula = props ~ Food + Temperature, data = Proportions, 
    subset = Instar == "5th")

Residuals:
      Min        1Q    Median        3Q       Max 
-0.047401 -0.027360 -0.002382  0.028821  0.046099 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)     0.514059   0.019202  26.771 1.22e-10 ***
FoodPoor diet  -0.079812   0.022172  -3.600  0.00485 ** 
Temperature30° -0.087631   0.022172  -3.952  0.00272 ** 
Temperature33°  0.007841   0.033259   0.236  0.81837    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.0384 on 10 degrees of freedom
Multiple R-squared:  0.7419,    Adjusted R-squared:  0.6644 
F-statistic: 9.579 on 3 and 10 DF,  p-value: 0.002749

The age structures of the larval populations are therefore significantly different between treatments, with a significant interaction between temperature and food quality. The 27 degree, good food populations have the most 5th instar larvae, and either increasing the temperature or changing to poor food decreases the relative abundance of 5th instar larvae in favour of younger animals. Food quality has little effect on age structure for the 30 degree populations, however, all of which have roughly similar age structures.

Time series analysis

Spectral analyses

#### Plotting spectral densities with bootstrap 95% limits derived from the 95th
#quantiles for spectral densities of 10,000 replicates of white noise drawn from
#a negative binomial distribution with the mean and shape parameter equal to 
#those for the time series in question

library(MASS)

adults<-read.csv("Dead adult counts.csv")

adults$temp<-factor(adults$temp)

rep<-ifelse(adults$box<7,1,ifelse(adults$box>12,3,2))

adults$food<-relevel(adults$food, ref="good")


adults<-data.frame(adults,rep)
rm(rep)

adults2<-subset(adults, week>10 & week<82)

attach(adults2)

sds<-tapply(total,box,function(x) sd(x, na.rm=TRUE))
means<-tapply(total,box,function(x) mean(x, na.rm=TRUE))

par(mfrow = c(3,4))
par(mar = c(1,1,1,1))

### Generate output data frame
Temperature <- c(rep(c(rep(c(27,30), each = 72),
                       rep(33, times = 36)), times = 2), 
                       rep(c(27,30), each = 72))
                   
Food <- c(rep(rep(c("Standard", "Poor", "Standard", "Poor", "Poor"), 
                  each = 36), times = 2), 
                  rep(c("Standard", "Poor", "Standard", "Poor"), each = 36))

Replicate <- c(rep(1, times = 180), rep(2, times = 180), rep(3, times = 144))
Freq <- numeric(0)
Spec <- numeric(0)
Quant <- numeric(0)

boxes<-c(1,2,3,4,6,7,8,9,10,12,13,14,15,16)



for(i in boxes) {

 
    
box_test<-i

shape<-theta.md(adults2$total[adults2$box==box_test], mu = means[box_test], 
                dfr = length(adults2$total[adults2$box==box_test]-1))

# temp2<-rnorm(71000, mean = means[2], sd = sds[2])
temp2<-rnbinom(710000, mu = means[box_test], size = shape)

sims<-matrix(data = ifelse(temp2>=0, temp2, 0), nrow = 10000, ncol = 71)

spec1<-function(x, smoother = 3) {
  
  smooth.spec = smoother
  temp<-spectrum(sqrt(abs(x)), span = smooth.spec, plot = FALSE)
  return(temp$spec)
  
}


spectra <- apply(sims, 1, spec1)

quant1 <- apply(spectra, 1, function(x) quantile(x, probs = c(0.95)))

smooth.spec<-3

temp1<-spectrum(sqrt(total[box==box_test & week <82]), span=smooth.spec, plot = FALSE)

Freq <- c(Freq, temp1$freq)
Spec <- c(Spec, temp1$spec)
Quant <- c(Quant, quant1)


# plot(temp1$freq, temp1$spec, type = "l")
# 
# points(temp1$freq, quant1, type ="l", lty = 2)
# 
# abline(v=1/6, lty=3, col="steelblue")

}

detach(adults2)

Spectral_output<-data.frame(Temperature, Food, Replicate, Freq, Spec, Quant)

Spectral_output$Temperature <- factor(Spectral_output$Temperature)
Spectral_output$Food <- factor(Spectral_output$Food)
Spectral_output$Food <- relevel(Spectral_output$Food, ref = "Standard")
Spectral_output$replicate <- factor(Spectral_output$Replicate)


par(mfrow = c(1,1))

temp_labels <-c("27" = "27°", "30" = "30°", "33" = "33°")
food_labels <-c("Standard" = "Standard food", "Poor" = "Poor food")
rep_labels <-c("1" = "Replicate 1", "2" = "Replicate 2", "3" = "Replicate 3")


p1 <- ggplot(data = Spectral_output, aes(x = Freq, 
                      y = Spec)) + geom_line(size = 0.4) + theme_bw()

p1 <- p1 + geom_vline(xintercept = 1/6, colour = "steelblue", 
                      linetype = "dotted", size = 0.5)

p1 <- p1 + geom_line(aes(x = Freq, y = Quant), linetype = "dashed", size = 0.4)

p1 <- p1 + xlab ("Frequency") + ylab("Spectral density")

p1 <- p1 + theme(panel.grid.major=element_blank(), panel.grid.minor=element_blank(), 
                 strip.background =element_rect(fill="white", colour = "grey80"))

p1<-p1+facet_grid(Replicate ~ Temperature + Food, 
                  labeller = labeller(Replicate = rep_labels, 
                  Temperature = temp_labels, Food = food_labels), 
                  scales = "free_x", space="free_x")
p1

Figure S6 Spectral densities for each population cage. Only those which persisted for the full period are shown. Dashed lines indicate approximate statistical significance levels: these show the 95th quantile for spectral densities of 10,000 replicates of white noise drawn from a negative binomial distribution with the mean and shape parameter equal to those for the time series in question. Dotted vertical lines indicate the frequency equivalent to a wavelength of six weeks. All populations at 27ºC and 30°C show strong peaks at a wavelength of six weeks, indicating regular cycles of one generation in length, apart from the 30°C populations with standard diet, which all show multiple, marginally significant peaks with little consistency. One 33°C population shows generation cycles but the other does not.

The data were square root transformed before analysis in order to normalise the error distributions. These smoothed spectral density plots show that in fact there are two important periodicities in the data - firstly a long-period cycle of between 35 and 25 weeks, as indicated by the peak on the left hand side of the graphs, and secondly short-period generation cycles of six weeks overlaid onto the longer term cycles. This is shown by the peak that in most cases aligns with the vertical dotted line, which indicates a frequency of 1/6 on each graph.

The locations of the peaks corresponding to the long-period cycles are given in table 2. Some populations show longer period cycles (approx. 35 weeks) and some shorter (approx 25 weeks). Note that because of the length of these cycles the spectral analysis will only fit peaks corresponding to periods of 72,36,24,18 etc. so the precise lengths of the long period cycles are hard to estimate.

popplot<-ggplot(data=adults2, aes(x=week, y=sqrt(total))) + labs(x="Time (weeks)",
                                    y="Square root number of dead adults" )

popplot<-popplot + geom_line(size = 0.25)

popplot<-popplot+facet_grid(rep~temp+food)+geom_smooth(method="loess",
                                    span=0.4, size=0.25, se=FALSE)

popplot<-popplot+theme_bw()+theme(panel.grid.major=element_blank(), 
                                  panel.grid.minor=element_blank())

popplot

Figure S7 The 18 population time series, square root transformed, with a fitted Loess smoother which makes the long-period cycles clear.

There is no apparent association between temperature and cycle length, but there is between food quality and cycle length. If we consider the 27 and 30 degree populations, all the poor food populations show short cycles (six out of six), whereas most of the good food populations show longer cycles (4 out of 6). This is not quite statistically significant:

cycles<-matrix(c(6,0,2,4), ncol=2)
rownames(cycles)<-c("Long.cycles","Short.cycles")
colnames(cycles)<-c("Poor.food","Good.food")
cycles

             Poor.food Good.food
Long.cycles          6         2
Short.cycles         0         4

fisher.test(cycles)

    Fisher's Exact Test for Count Data

data:  cycles
p-value = 0.06061
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
 0.9097849       Inf
sample estimates:
odds ratio 
       Inf 

Regarding the generation cycles, the spectral densities reveal strong peaks corresponding to a period of 6 weeks in all of the 27 degree populations - these populations have dynamics that are essentially dominated by a long period cycle with the shorter period generation cycles overlaid onto them, and there is little indication of other important periodicity. For the 30 degree populations, there are strong peaks at six weeks in all three of the “poor food” populations although the picture is not as simple for these populations and population 14 has a notable third peak corresponding to a period of roughly 10 weeks. The 30 degree populations given “good” food have complex spectra with multiple peaks and little indication of strong generation cycles.

Of the 33 degree populations which persisted, one (population 6) had little to indicate generation-length cycles, whereas the other (population 12) had a strong peak at a frequency corresponding to a six-week period.

Wavelets

Wavelet analysis (Cazelles et al. 2008) allows us to examine whether any periodic behaviour exhibited by our populations changed during the 82 weeks they were maintained, and is used here to try to understand the differences in periodic behaviour between the 27ºC and 30ºC populations. Because the long-period cycles in these populations have wavelengths which are greater than 25% of the total length of the time series, these were removed from the data by smoothing with a loess function of span 0.4 prior to analysis. Wavelet transformations were carried out using the WaveletComp package (Roesch & Schmidbauer 2014) using a Morlet wavelet and plotted with 5% significance contours calculated using the default (white noise) options.

library(WaveletComp)

######## Plotting the wavelets without the population data.




par(mfrow = c(3,5))

par(mar=c(1,1,2,1)) #Margin sizes

to.plot<-c(1,2,3,4,6,7,8,9,10,12,13,14,15,16) 
#Only plot the 27 and 30 degree populations plus 2 33 degree ones which persisted

for (i in to.plot) {
    
    box<-adults2$total[adults2$box==i]
    weeks<-adults2$week[adults2$box==i]
    temp<-loess(box~weeks, span=0.4)
    #plot(temp$residuals, type="l", main=paste("Population ", i), xaxs="i")    
    #Plot the residuals after a loess fit - 
    #in other words the data with the long wavelength cycles removed


    mydata<-data.frame(x=box)
    my.w2 = analyze.wavelet(mydata, "x",          
                            #Calculate wavelet power spectrum 
                            loess.span = 0.4,   
                            #Remove long wavelength cycles before analysis
                            dt = 1, dj = 1/250,
                            lowerPeriod = 2,    
                            #Check for periods between 2 and 24 weeks
                            upperPeriod = 24,
                            make.pval = T, n.sim = 10)     
                            #Use more simulations for the significance test 
                            #than the default which gives smoother contours
                            
    wt.image(my.w2, color.key = "interval",     
             #Plot wavelets on an interval scale
                    n.levels = 250,color.palette = "rainbow(n.levels, start=0, end=0.7)",
                    legend.params = list(lab = "wavelet power levels", mar = 4.7), 
                    graphics.reset=FALSE,    
             #Need this so that it doesn't just overplot in the same place
                    plot.ridge=FALSE,      
             #Plot "Ridge of power"
                    lvl=0.8,              
             #Only plot ridge for reasonably high powers 
                    siglvl=0.01,         
             #Set p-value for significance contour to 0.05 rather than the default 0.1
                    label.time.axis = FALSE, label.period.axis = FALSE,
                    plot.legend=FALSE)     
    
    mtext(text = paste("Temp = ", adults2$temp[adults2$box == i][1], 
                       " Food = ", adults2$food[adults2$box == i][1]), cex = 1.3)  
            #Add text above plot giveing temperature and food quality
    }


    
    par(mfrow=c(1,1))

Figure S8. Wavelet power spectra of population time series. The x-axis gives time in weeks and the y-axis indicates the period. The power is indicated by colour, with blue indicating low power whereas green, yellow and red show increasing power, and the paler regions at the left and right hand sides of the plots show the regions where effects associated with the beginning and end of the time series reduce certainty. The time series identified by spectral analysis (Fig. 4) as showing 6-week generation cycles, including all of the 27°C populations, the 30°C poor diet populations and the second replicate of the 33°C poor diet treatment all show areas of high power at a period of 6 weeks, with those with particularly strong 6-week peaks in the spectral analysis showing ‘ridges’ of power across the time series, indicating that generation cycling is occurring for most of the period that the populations were maintained. The 30°C populations maintained on standard diet, by contrast, show brief periods of regularity at a variety of periods, and only the second one shows any peak of power near the 6-week periodicity found in the other time series. The 33°C population which does not show generation cycles (the first one) has a very brief period of cycling at the start of the time series after which there is no real periodicity. Populations which became extinct during the experiment are not shown.